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•-^ In most models of structure formation, the primordial den- 
sity fluctuations are assumed to have Gaussian statistics, 
statistical properties of the density field in this case 
^ are then completely specified by the two-point correlation 
\—^ function or equivalently by its Fourier transform, the power 
^ spectrum. Our knowledge of these quantities has improved 
greatly within the past few years with the analysis of cluster- 
ing in several large, recently completed galaxy surveys (e.g. 
Davis & Peebles 1983, Maddox et al. 1990, Efstathiou et al. 
1990, Saunders et al. 1991, Strauss et al. 1992, Loveday et 
al. 1992, Fisher et al. 1992, Vogeley et al. 1992, Baugh & 
Efstathiou 1993, 1994a; Fry & Gaztafiaga 1994) and with 
the detection of anisotropics in the microwave background 
radiation (e.g. Smoot et al. 1992, Hancock et al. 1994). How- 
ever, we also need to measure the higher order moments of 
the density field to test for consistency with a Gaussian dis- 
tribution. 

For a Gaussian density field, all the connected moments 
of the distribution function of the fluctuations are zero. As 
the fluctuations grow under the influence of gravity, the vari- 
ance of the density field increases; as the variance approaches 
unity, the distribution of density fluctuations becomes asym- 
metrical, developing non-zero higher order correlations, i.e. 
skewness, kurtousis and so on. Non-zero skewness has been 
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detected in several galaxy catalogues (e.g. Groth & Peebles 
1977, Saunders et al. 1991, Bouchet et al. 1993, Gaztafiaga 
1992, 1994). It is possible however, in cosmological mod- 
els that contain primordial structures that are nonlinear, 
such as cosmic textures or strings, that the initial density 
field is non-Gaussian (Silk & Juskiewicz 1991, Weinberg & 
Cole 1992). The question that then needs to be answered 
is whether the observed skewness can be explained by the 
gravitational evolution of an initially Gaussian density field, 
or whether some degree of primordial skewness is required. 

In order to answer this question, we need to follow the 
evolution of the density field into the nonlinear regime, de- 
fined by the density contrast becoming on the order of unity, 
8p/p ~ 1, and larger. Rather than trying to estimate the un- 
derlying distribution of densities directly, we shall measure 
the moments of the density field using the statistics of counts 
in cells (see for example Peebles 1980, §36). The method con- 
sists of dividing up the volume of space under consideration 
into cells of side R and calculating the moments of the num- 
ber counts of objects in these cells. The connected moments 
provide then a volume-averaged measure of the N-point cor- 
relation functions, fj^. 

Analytically, Peebles (1980) used second-order pertur- 
bation theory (PT) to calculate the skewness at a point in 
the density field obtaining the simple result, ^3 = ^3/^2 — 
34/7. In practice, the skewness is measured averaged over 
some finite volume and one has to study the multipoint 
contribution (Fry 1984), and then smooth it (Goroff et al. 
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1986). For a finite splierical cell, Juskiewicz et al. (1993) 
computed the the skewness for a scale free initial power spec- 
trum P{k) ~ and found a dependence upon the spectral 
index n, S3 = 34/7 — (ra + 3). Frieman & Gaztafiaga (1994) 
estimated ^3 numerically as a function of scale R for a CDM 
spectrum. Recently, these results have been generalized for 
an arbitrary power spectrum and higher order moments by 
Bernardeau (1994a, 1994b). 

Baugh & Efstathiou (1994b, hereafter Paper 1), com- 
pared the nonlinear evolution of the power spectrum of den- 
sity fluctuations modelled by A'-body simulations with the 
predictions of second-order perturbation theory to estimate 
the range of validity of the perturbation theory results. In 
this paper we shall extend this comparison to higher order 
moments of the density distribution. 

There have been several studies of skewness and higher 
order correlations in N-body simulations, dealing both with 
scale free and CDM models (e.g. Efstathiou et al. 1988, 
Weinberg & Cole 1992, Bouchet & Hernquist 1992, Lahav et 
al. 1993, Fry, Melott & Shandarin 1993, Colombi et al. 1994, 
Lucchin et al. 1994, Bernardeau 1994a). In this paper, we 
make a much more detailed comparison with perturbation 
theory than attempted in previous studies. Different sized 
CDM simulations both in terms of the number of particles 
(from N = 32768 up to TV = 10^) and the computational box 
length (up to Lb = 300 Mpc) are used to check finite 
volume and shot noise effects. Most of the previous analyses 
have used simulation boxes of side Lb < 50 Mpc which, 
as we show below, are subject to very large systematic fluc- 
tuations in Sj. Here, each simulation is represented by an 
ensemble of 10 different realizations of the random phases 
which makes it possible to estimate precise and realistic er- 
rors, while in previous analysis very crude errors or no errors 
at all have been given. In each ensemble we follow all the 
higher moments from 2nd to 10th order for a wide range 
of length scales (up to 75 Mpc) as they evolve in time 
(for up to 12 expansion factors), whilst previous analyses 
have focused on smaller scales, lower statistics, J < 4, and 
just one output time. Moreover, a detailed comparison of 
the amplitudes Sj for J > 3 with PT was not possible until 
just very recently (Bernardeau 1994b). 

The layout of the paper is as follows. In Section 2 we give 
a brief description of the counts in cells method and show 
how well we can reproduce the variance of the density field 
using this technique. We give a systematic study of the mea- 
surement of the moments of the density field from numerical 
simulations in Section 3. In addition, we examine the vari- 
ous schemes for estimating the errors on the moments that 
have been used in the literature and compare these with the 
ensemble errors. The comparison between the moments of 
the density field measured from the simulations and the pre- 
dictions of perturbation theory is made in Section 4. Finally, 
we summarise our conclusions in Section 5. 



2 THE COUNTS IN CELLS METHOD 

The calculation of the moments of a particle distribution 
and their relation to the corresponding moments of the con- 
tinuous density field have been discussed at length in the 
literature (see for example Peebles 1980, §36 or Gaztafiaga 



1994). We shall just quote here some of the formulae for the 
moments that we use in this paper. 

If the volume containing the particle distribution is di- 
vided up into cells of comoving radius R, the j"* central 
moment of the counts of the mass distribution is given by 

M 

mj(R) = J2iN'-W, (1) 

8 = 1 

where N, is the number of particles in the i*^ cell, N is the 
mean number of counts for cells of size R and the summation 
is over the M cells in the volume. 

Using the notation of Gaztafiaga (1994), the volume 
averaged connected correlation functions, fij, can be written 
in terms of the mj 

fi2 = m2 
II3 = ms 

fii = nii — 3m2. (2) 

The discreteness of the particles gives rise to a extra 
contribution to the moments of equation 2, which becomes 
significant on scales for which the number density of the 
particles is around unity. If the particles have been drawn 
at random from some underlying parent distribution, the 
noise effects can be corrected for using the Poisson shot 
noise model. Applying this correction leads to the follow- 
ing expressions for the connected moments up to J = 4: 

k2 = m2 — N 

ks = ms — 3m2 -|- 2N 

ki = nii — 3m2 — 6m3 -|- llm2 — 6N. (3) 

The volume averaged J— point correlation function, 
S,j{R), is defined as 

Jj = -L. [ dn ...drjW(ri) ... W(rj)ij(ri,...,rj), (4) 
J 

where the density fluctuations have been smoothed over a 
window function FF(x) with volume Vw. These moments can 
be obtained from equations (2) and (3) simply by dividing by 
N'' ; for example the volume averaged two-point correlation 
function can be written as 

?2(-^) = m2/N (uncorrected) (5) 

^^{R) = m2/N -1/N (shot noise corrected). (6) 

We shall employ the notation throughout the rest of the 
paper, indicating in each case whether or not any discrete- 
ness corrections have been applied. 

2.1 A test of the counts in cells method 

We first compare the measurement of the second moment 
using the counts in cells technique, with an estimate of the 
variance obtained using the power spectrum of the density 
fluctuations. To make the best possible measurement of the 
second moment we use the largest simulations that are avail- 
able to us, i.e. the 100'^ particle simulations in a box of 
300/t"^Mpc ( see §3 and Table 1). We use spherical cells to 
obtain an estimate of the second moment for each simulation 
in ensemble B. We then average over the 10 members of the 
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Figure 1. The volume averaged 2-point correlation function ^2 
for different output times in the 100'^ particle simulations as a 
function of the comoving radius of the sphere. Squares with error 
bars show the estimation of the variance using the counts- in- cells 
technique, whilst the dashed lines correspond to the Fourier trans- 
form of the estimated power spectrum P{k). The solid line shows 
the variance computed from equation (7) when P{k) is extrapo- 
lated to small and large k as described in the text. 



ensemble of initial random phases to find the la variance 
on the mean. The second moment measured in this way is 
shown by the symbols in Figure 1. 

Using the power spectra measured in Paper I for en- 
semble B, we can plot an another estimate of the variance 

^2 = ^ fj Akk^P{k)W\kR), (7) 

where W(kR) is the Fourier transform of the spherical win- 
dow with radius R 
3 



W{kR) 



(kR)' 



[sm(kR) - kRcos(kR)]. 



(8) 



The simulation box and the number of particles used set 
upper and lower limits, ki and k2 respectively on the range 
of wavenumbers for which we can represent the theoretical 
power spectrum of density fluctuations (c/ Appendix): 



ki = 27r/i 



k2 = 27r/i 



(9) 



where L is the length of the simulation box and Npar is the 
total number of particles. The second moment computed 
from equation (7) with these limits is shown by the dashed 
lines in Figure f. 

There is excellent agreement between the two estima- 
tions of the second moment down to scales of ~ 4/t~^Mpc. 
For smaller scales the values of {2 from equation (7) are lower 



than those measured using counts in cells. This is caused 
by the truncation of P{k) beyond the Nyquist frequency; 
P{k) = for k > k2 (dashed line). If we extend the mea- 
sured power spectrum beyond the Nyquist frequency with 
a linear extrapolation of the last few points we obtain the 
curves shown by the solid lines in Figure 1. These curves 
show a better agreement with the counts in cells results on 
scales less than 4/t~"'Mpc. 



3 COUNTS IN CELLS ANALYSIS OF THE 
SIMULATIONS 

In this Section we present a systematic examination of the 
practical considerations involved in performing a counts in 
cells analysis of simulation data. We address the effects that 
the following have upon the calculation of the moments: 

i) The intial arrangement of particles that is perturbed to 
set up the initial density fluctuations. 

ii) The discreteness corrections and the finite volume of 
the simulation box. 

iii) The shape of the smoothing window; spherical or cubi- 
cal. 

iv) Different schemes for estimating the errors on the mo- 
ments. 

In order to investigate these points, we have run several 
new ensembles of A'-body simulations to use alongside those 
of Paper 1. All the models used are of a standard Cold Dark 
Matter (CDM) universe (F = Q/t = 0.5) evolved using a 
P"^ M code. The initial linear power spectrum of the density 
field is that of Bond & Efstathiou (f 984) : 



P{k) oc 

where 

J/ 

a 
b 

c 



[1 + (ak + (hkfl^ + (ckfYfl^ 



f.f3 



(10) 



6.4/(Q/j^)Mpc 
3.0/(Q/t^)Mpc 
f.7/(Q/j^)Mpc. 



This form for the power spectrum applies for scale-invariant 
CDM universes that have low baryon densities, Q_b ~ 0.03. 
The parameters of these simulations are listed in Table f; 
the simulations used in Paper 1 are Ensembles A and B. 
The final column gives the initial particle arrangement from 
which the particles are displaced to set up the theoretical 
spectrum of density perturbations. The glass initial condi- 
tions refer to a particle distribution that is subrandom with 
no regular structure, with the particles avoiding one another. 
A description of this distribution and its production is given 
in the Appendix, along with a power spectrum comparison 
of the clustering of mass in the ensembles containing 32'^ 
particles. We shall refer to the ensembles of simulations by 
the number of particles that they contain and the pattern 
of particles that was displaced to set up the initial density 
fluctuations. 

We also list the length scales corresponding to the 
Nyquist frequency of the particle distribution (Nyquist 
length) and that for which the number density of particles 
in the simulations is unity (shot length). The Nyquist length 
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represents the shortest length scale on which the theoretical 
power spectrum of density fluctuations can be represented 
in the initial conditions of the simulation (see Appendix for 
a discussion of this point); in the grid and glass ensembles 
this length scale is given by the Nyquist frequency of the 
particle grid 

kNyg = (27r/i) (11) 

For the random ensemble (D), the theoretical spectrum 
is truncated on even larger scales (see the Appendix). The 
shot length gives the scales for which the discreteness of 
the particles begins to have a significant effect upon the 
moments measured in the simulations. 

We estimate the moments of the particle distribution in 
the ensembles by dividing the simulation box up into cubical 
cells. The smallest scale that we examine in the simulations 
is half the shot length; the largest scale corresponds to cubi- 
cal cells of side equal to one half the length of the simulation 
box. The increment in volume of the cells is by a factor of 
~ 1.5, so that different parts of the power spectrum of the 
density fluctuations are sampled by the cell window func- 
tion. We also measure the counts in spherical cells, which 
are also displaced relative to the cubical cells in order to 
fully sample the mass distribution. 

We examine the density field modelled by the simula- 
tions at different epochs. We will quantify the evolution of 
the density field in terms of the linear variance measured 
in spheres of radius 8/t~^Mpc, denoted ag, calculated using 
equation 7, with limits ki = and k2 = oo. Occasionally we 
shall refer to the epoch using a redshift, which relies upon 
identifying a reference epoch in the simulations, usually that 
for which erg = 1, as the present day. 

3.1 Initial particle arrangement 

We compare the variance measured in the simulations 
with 64'^ particles started from grid and glass initial arrange- 
ments of particles (ensembles A and F respectively) in Fig- 
ure 2. The variances shown are the average over 10 simula- 
tions in each ensemble and are plotted with the la errors 
on the mean. We show the variance at two epochs in the 
evolution of the density field, corresponding to ag = 0.50 
and ag = 1.00 . At the earlier epoch in the evolution of 
the density field, there is some discrepancy in the variances 
measured in the two ensembles, with the simulations started 
from a glass initial pattern showing a higher variance. At 
later times, and for larger scales there are negligible differ- 
ences between the variances measured from the two ensem- 
bles. 

Figure 2 illustrates the difficulty in assessing whether 
the deviation of the measured variance from the linear the- 
ory prediction is the result of genuine non linear evolution 
of the density or is merely due to discreteness noise. 

3.2 Discreteness and finite volume effects 

In running a A'-body simulation, we are trying to ap- 
proximate a continuous density field with a discrete set of 
masses contained inside a box of finite size. The limitations 
of this approximation introduce artificialities into the mod- 
elled density field on both large and small scales. 
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Figure 2. A comparison of the variance measured in the simu- 
lations of ensembles A (64'^ particles started from a grid initial 
pattern - filled squares) and F (64'^ particles started from a glass 
initial pattern - open squares). The output times shown corre- 
spond to (Tg = 0.25 and ag = 1.00. The solid lines show the 
variance predicted by linear perturbation theory at these epochs. 



3.2.1 Large scale effects 

Due to the finite size of the simulation box, the power spec- 
trum of the density fluctuations is truncated artificially on 
scales larger than the box. Also, as the simulation evolves, 
the lowest order Fourier modes may eventually become non- 
linear. When this stage is reached, the simulation has to be 
stopped, because the nonlinear interactions can no longer be 
modelled accurately on length scales equal to the size of the 
box (Davies et al. 1985, Paper 1). 

We examine finite volume effects by plotting the vari- 
ance measured from the simulations of ensembles A, B, and 
C on the same axes in Figure 3. Ensembles A and C have 
box sizes of 180/t~^Mpc, whilst the simulations of ensemble 
B are in boxes of side 300/t~^Mpc. There are only small dif- 
ferences at the largest scales, R < LbI^, plotted for each 
simulation, indicating that the effect is small. 

To simulate the finite volume effect, we have estimated 
^2 by truncating the linear P{k) at ki = 21: 1 Lb, which 
corresponds to the largest scale in each ensemble to mimic 
the lack of power at larger scales. We find that the effects 
are not very important at i? < LbI^ (c/ figure 9). Thus by 
using R = Lb /8 for the largest cell radius our results are 
unaffected by the the finite volume of the simulation. 



3.2.2 Small scale effects 

On scales much smaller than the size of the simulation box, 
the discreteness of the particles becomes important. An ex- 
tra variance or shot noise is present on scales around the 
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Table 1. #-body simulation parameters 



Ensemble 


Number of 


number of 


mesh size 


box size 


Nyquist length 


shot length 


intial particle 




simulations 


particles 










arrangement 


A 


10 


64^ 


128^ 


180 


5.6 


2.8 


grid 


B 


10 


100^ 


256^ 


300 


6.0 


3.0 


grid 


C 


10 


32^ 


64^ 


180 


11.3 


5.6 


grid 


D 


10 


32^ 


64^ 


180 


30.0 


5.6 


random 


E 


10 


32^ 


64^ 


180 


11.3 


5.6 


glass 


F 


10 


64^ 


128^ 


180 


5.6 


2.8 


glass 
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Figure 3. A comparison of the variance measured in the sim- 
ulations of ensembles A (64'^ particles - open squares), B (100'^ 
particles - open stars) and C (32'^ particles - filled squares), all 
started from a grid pattern, at output times corresponding to 
(Tg = 0.50 and (Tg = 1.00. The solid lines show the variance pre- 
dicted by linear perturbation theory at these epochs. 



mean interparticle separation because the approximation 
that the density field is continuous breaks down on these 
scales. In addition to this, the theoretical power spectrum is 
truncated in the initial conditions at the Nyquist frequency 
of the particle grid (see Appendix). Both these factors affect 
the accuracy with which the small scale density fluctuations 
can be represented. The glass initial arrangement shows no 
significant increase in dynamic range in length over the grid 
simulations, for the relatively flat power spectrum of the 
CDM model. However, a glass initial pattern is more useful 
in the modelling of fluctuations with a steeper spectrum, in 
which voids in the particles distribution are more prominent, 
as is the case in the Hot Dark Matter scenario (White 1993, 
private communication). 

Figure 3 compares the discreteness effects on small 
scales in simulations with different numbers of particles. The 



simulations in ensembles A and B have approximately the 
same Nyquist frequency, corresponding to a length scale of 
~ 6/t"^Mpc. The shot noise in the simulations of ensemble 
B is however a factor of 4 lower than that present in the sim- 
ulations of ensemble A. The simulations with 32'^ and 64'^ 
particles show variances that are in very good agreement 
with that measured in the 100'^ particle simulations, down 
to the length scale corresponding to their respective Nyquist 
frequencies. On smaller scales than this, the shot noise dom- 
inates and leads to discrepancies between the measured vari- 
ances, particularly at earlier epochs in the evolution of the 
density field. We defer a discussion of whether or not it is 
possible to model the shot noise present in the simulations 
to the end of this Section. 



3.3 The shape of the smoothing window 

In Figure 4, we plot the ratio of the variance measured us- 
ing cubical cells to that measured with spherical cells at 
two different output times in the 100'^ particle grid simula- 
tions. This is the variance without any correction for particle 
discreteness, corresponding to equation (5). The error bars 
show the 1(7 variance on the mean averaged over the simu- 
lations of the grid ensemble. 

To compare the measurements made with different win- 
dow shapes, we have rescaled the dimensions of the cubical 
cells I to the radius of a sphere that would have the same 
volume. Re = (3/47r)^/^l. There are significant differences 
in the moments measured using different cell geometries for 
the initial conditions in the simulations. This is due to in- 
terference between the cubical cells and the initial grid po- 
sitions of the particles, particularly on scales for which the 
shot noise is important. At later times in the simulations, 
the difference between the measurements becomes smaller, 
but is still significant. The choice Re = (3/(47r))"'^'^l, though 
natural, is not necessarily the right one. For a power-law 
power spectrum it is easy to find the value of R that corre- 
sponds to I by just estimating { numerically. If the size of 
the cubical cell is scaled to i? ~ 1.025i?e we find excellent 
agreement, i.e. top of Figure 4. 

To ensure full coverage of the simulation box, the spher- 
ical windows are laid down so that they slightly oversample 
the point distribution, while the cubical cells exactly sample 
the box just once. Thus, on large scales one has a relatively 
small number of cubical cells, which explains why the ratio 
in Figure 4 shows larger error bars for R > 20/t~^Mpc. 

Thus, we shall use spherical cells henceforth in this pa- 
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Figure 4. The ratio of the variance measured using cubical cells 
to that measured with spherical cells for the large simulations 
(Ensemble B) at different output times. In the top plot we have 
multiplied the radius of the equivalent spherical cell by 2.5%. 



per, both to avoid possible spurious effects due to interac- 
tions with the particle grid and because the perturbation 
theory calculations are simpler for spherical windows. 



3.4 Error estimation 

The initial conditions for A'-body simulations require the 
Fourier transform of the density field to be set up with an 
amplitude and a random phase at each grid point in Fourier 
space. Rather than using just one simulation, which could 
be subject to spurious random fluctuations, we need to aver- 
age over several realisations of the density field before com- 
paring the predictions of the simulations with observation 
or analytic calculations (Efstathiou et al. 1985). For practi- 
cally all purposes, averaging over 10 simulations appears to 
be sufficient, as illustrated by Figure A3 in the Appendix, 
which shows how small the la errors in the power spectrum 
of density fluctuations are averaged over ensembles of this 
size. 

To calculate the variance on the second moment for- 
mally, a knowledge of the fourth moment is necessary 
(Kendall and Stuart 1977) 



Var({2 



2{3 + £.2 



M 



(12) 



where M is the number of cells in which the counts are mea- 
sured. The best way however, to estimate the uncertainty in 
the variance is to average over the simulations in the ensem- 
ble and to calculate the variance on this mean. The error 
bars that we show in this paper are the la variance on the 
mean J"^— moment computed in this way. 

In Figure 5 we compare the values of {2 for each indi- 
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Figure 5. Ratios of the correlations ^2 in each of 10 individual 
simulation to the average in the ensemble A (upper graph) and 
ensemble B. The errorbars represent the 1-sigma interval. 



vidual simulation. For clarity, we have scaled {2 to the mean 
in the ensemble. The errors in the larger simulation are sig- 
nificantly smaller. It is interesting to note that the values 
of £,2{R) in a given simulation do not necessarily fluctuate 
around the mean function < £,2{R) > but can be signifi- 
cantly shifted at all scales with respect to the mean by up 
to ~ 5% in the Lb = 300 Mpc ensemble and ~ 10% 
in the Lb = 180 Mpc one. These overall shifts are pro- 
duced by density fluctuations on scales comparable to the 
box size which are not properly represented in each individ- 
ual simulation. 

In Figure 6 we compare the values of ^3 = S.s/S.2 for 
each individual simulation within the Ensembles. Again, the 
errors in the larger simulation are significantly smaller and 
the value of S3 measured in a given simulation does not nec- 
essarily fluctuate around the mean < S3{R) >, but can be 
significantly shifted at all scales with respect to the mean. By 
using just one simulation and underestimating the sampling 
error, there is a large probability of missing the perturbation 
theory agreement we find below (see Figure 11 and Figure 
13). We have also estimated the variance of S3 in smaller vol- 
umes. The fluctuations are typically 2-3 times larger for a 
sample of Lb ~ 100 Mpc than in the Lb = 180 Mpc 
one. The typical sampling error in a Lb = 100 Mpc sam- 
ple £at ~ 10-20 h-^ Mpc) is 10% for J^, 25% for {3 and 35% 
for These factors increase even more for smaller samples 
and makes the Sj estimations very unreliable in small sim- 
ulations or small catalogues, i.e. for typical redshift samples 
(see also Gaztafiaga 1994). 

It is instructive to compare the ensemble errors that we 
employ in this paper with the various schemes for error es- 
timation that have been used in the literature, when only 
one simulation has been analysed. This will allow a compari- 
son of the results presented in this work with those reported 
elsewhere and more importantly, it will allow the relative 
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Figure 6. Hierarchical skewness Ss in each of 10 individual sim- 
ulations in Ensemble A (upper graph) and Ensemble B. The er- 
rorbars are obtained from the 1-sigma errors on the mean values 
ofig and 

merits of these schemes when applied to real data sets to be 
evaluated. 

We shall consider four error estimation schemes; (i) 
splitting the simulation up into zones and averaging over 
the zones (see for example Maddox et al. 1990 for an ap- 
plication to the angular correlation function of galaxies), 

(ii) bootstrap resampling (as used by Lahav et al. 1991; see 
Ling, Frenk & Barrow 1986 for details of implementation) , 

(iii) regridding of the counts and recalculation of the second 
moment, and finally (iv) averaging over random subsets of 
cells for a given cell size (Gaztafiaga 1994). 

Splitting the survey into zones and using the variance 
between the moments measured for the cells in each zone 
provides a conservative estimate of the errors because the 
cosmic variance is larger for the zones than it actually is for 
the full simulation. 

Bootstrap resampling involves making new samples 
from the original data and measuring the fluctuations over 
these samples. The new sample is made by drawing at ran- 
dom from the list of mass points in the simulation with re- 
placement (i.e each point can be chosen more than once) un- 
til the bootstrap sample contains the same number of points 
as the simulation. Hence at the positions occupied by mass 
particles in the simulation, there will be or 1 particles 
in the bootstrap sample, and with decreasing probability 
2, 3, 4, etc., particles. 

By regridding the simulation and redoing the counts in 
cells analysis, we hope to average over the cases where a 
dense cluster of mass points falls entirely within one cell or 
straddles the boundary between two cells. Such events could 
significantly alter higher moments and hence bias a formal 
calculation of the error on the variance. 

Scheme (iv), averaging over random subsamples of cells 
drawn at random from the simulation corresponds to a nu- 
merical estimation of the formal variance, i.e. equation (12), 



Table 2. A comparison of different schemes for estimating the 
errors on the measured variance. The first number in each column 
is the measured variance. The figure in brackets is the percentage 
error on the variance found using each method. 



method 


length scale [h "'Mpc) 




inn 


40. U 


ensemble 


0.558 ± 1.6% 


0.0137 ± 16% 




0.573 ± 7.5% 


0.0108 ±45% 


(i)(b) 


0.550 ± 7.1% 




(ii)(a) 


0.561 ± 0.2% 


0.0170 ± 0.6% 


(n)(b) 


0.561 ± 0.1% 


0.0170 ± 0.1% 


(iii) (a) 


0.570 ± 3.2% 


0.0117 ± 25% 


(iii)(b) 


0.563 ± 1.4% 


0.0169 ± 11% 


(iv)(a) 


0.562 ± 6.4% 




(iv)(b) 


0.557 ± 7.7% 





and has a number of attractive features. Firstly, the prob- 
lem of a larger cosmic variance when a simulation is split up 
into smaller volumes is avoided because the cells are drawn 
from the entire simulation volume. Also, the method does 
not require periodic boundary conditions to be implemented 
in its fullest sense as with regridding. 

In all cases, we define the variance on the mean of the 
second moment, ({2) by 

I 

i = l 

where / is the number of trials or measurements of the sec- 
ond moment, e.g. the number of zones used or the number 
of bootstrap resamplings made. 

We present a comparison of these error estimation 
schemes in Table 2, using the final output time of a simu- 
lation from the grid ensemble with 32'^ particles, the size of 
simulation that is most commonly analysed in the literature. 
Two scales are examined; the length scale corresponding to 
the Nyquist frequency for ensemble C, i? ~ 10/t~^Mpc and 
the largest scale cells that we use, R = 45/t~"'Mpc. The 
former scale is used because this is the scale on which non- 
linear effects start to be dominated by discreteness effects. 
The latter scale was chosen because the counts in cells of 
this size will show large fluctuations due to the small num- 
ber of independent cells in the simulation volume. We show 
the mean moment and the percentage error on the mean for 
each method. 

The details of the various error estimation schemes pre- 
sented in Table 2 are as follows: 

i) The simulation was split up into (a) four zones and (b) 
eight zones; the table entry gives the mean and variance 
on the second moment computed from the cells in each 
zone. 

ii) (a) 100 and (b) 1000 bootstrap resamplings of the sim- 
ulation mass points were made, with replacement. The 
value for {2 given is calculated for the real simulation; 
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Figure 7. (a)The variance measured using spherical cells for the simulations of the grid (filled squares) , random (stars) and glass (unfilled 
squares) ensembles. The output times plotted correspond to (Tg = 0.24, 0.51 and 1.0. The linear perturbation theory variances are shown 
by the solid lines, (b) The variances of (a) corrected for the shot noise of the initial unperturbed particle distribution, (c) The variances 
of (a) corrected using the evolved noise variance as described in the text. 



the error is the variance over the estimates of {2 ob- 
tained from the resamplings 

iii) Regridding of counts using the periodic boundary con- 
ditions of the simulation; (a) using 4 regriddings, (b) 
using 8 regriddings. 

iv) Random subsets of cells drawn from the whole simula- 
tion volume; (a) 4 random subsets, (b) 8 random sub- 
sets. 

Table 2 shows that splitting the simulation up into zones 
gives the largest errors, as expected. The bootstrap tech- 
nique gives errors that are unrealistic. Simply increasing the 
number of bootstrap samples reduces the magnitude of the 
error. As there is no constraint on the number of bootstrap 
samples that one should make, it is difficult to see how this 
method can give us a reliable estimate of the error. Regrid- 
ding techniques give the smallest variances on the mean. 
Regridding is easier to motivate than making bootstrap sam- 
ples, because it is possible to have a scenario in which a dense 
clump of points straddles the boundary between two cells 
and is hence underweighted, or vice-versa is overweighted if 
it lies entirely within one cell. 

A further disadvantage of the above methods that in- 
volve examining the scatter between subgroups of the data 
is that these methods are not applicable for the largest cells 
used. On these scales, as mentioned above, the density fluc- 
tuations in the simulation box are poorly represented. The 
only way to properly measure the fluctuations on scales ap- 
proaching the box length is to average over an ensemble of 
simulations. 

3.5 Models for the shot noise 

Figure 3 shows that on small scales, particularly at earlier 
times in the evolution of the density field, there are signif- 
icant discrepancies between the variance measured in sim- 



ulations that contain different numbers of particles. These 
scales are also the scales at which the measured variances be- 
gins to move away from the prediction of linear perturbation 
theory. If we are to make a comparison of the predictions of 
perturbation theory with the results of the simulations down 
to these length scales, we need to assess how much of this 
deviation at scales less than ~ 5/t~^Mpc is due to true non- 
linear evolution and how much is the result of shot noise. 

In Figure 7(a) we plot the second moments measured 
in the 32'^ particle simulations started from grid, glass and 
random initial patterns of particles. No correction has been 
made to the second moment obtained from the ensembles. 
The moments of the grid and glass ensemble deviate away 
from the linear theory prediction on scales < 10/t~^Mpc in 
the initial conditions. The second moment of the random 
ensemble sits above the linear theory prediction; the initial 
pattern of particles has significant shot noise on all scales 
for this ensemble. 

The Poisson shot noise model is only strictly applicable 
when the points under consideration have been selected at 
random from a parent population; this is not a valid cor- 
rection to use for the moments of the dark matter particles 
in a simulation. It is a better approximation however for 
galaxies, if they trace the underlying dark matter distribu- 
tion randomly, though not if the bias between the dark and 
luminous matter is more complicated, for example if it is a 
function of density threshold (Kaiser 1984), non-linear (Fry 
& Gaztafiaga 1993) and/or scale dependent (see Bower et al. 
1992, Frieman & Gaztafiaga 1994 and Gaztafiaga & Frieman 
1994). 

As the displacements of the particles made to set up the 
initial density fluctuations are small, we can approximate 
the total variance {2 by the direct sum of the fluctuations 

—SN 

arising from the unperturbed arrangement of particles {2 > 

J (J 

with the fluctuations imposed by the displacements {2 • 
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-SN 



(14) 



We subtract {2 measured lor each initial pattern ol 
particles Irom the corresponding second moments at the se- 
lected output times in Figure 7(b). The agreement with lin- 
ear theory is now greatly improved, at least down to the 
scale corresponding to the Nyquist Irequency lor the glass 
and grid ensembles. Hence, at least in the initial conditions 
when the displacements ol the particles are small, we can 
correct lor the shot noise present using equation(14) (c/ 
also §4.1 below). However, there are still discrepancies at 
later times between the corrected moments lor the ensemble 
started Irom random initial positions and those started Irom 
grid and glass arrangements. As the simulation evolves the 
clustering ol the points grows. This evolution affects both 
the shot noise and the initial displacements. 

We attempt to model the evolution ol the shot noise 
in the simulations by doing a separate N-body simulation 
which only contains the shot noise contribution. We impose 
a white noise power spectrum onto the initial particle ar- 
rangement allowing the resulting perturbations to grow un- 
der gravity. The white noise spectrum gives the particles 
the same root mean squared displacements away Irom their 
initial positions that they exhibit in the CDM simulations, 
without introducing any coherent structures. One simulation 
ol this noise was run lor each set ol initial unperturbed posi- 
tions, grid, random and glass. We measured the variance in 
the point distributions at expansion lactors matching those 
used in ensembles C-E, using spherical windows. The errors 
are estimated by regridding the point distributions 10 times. 
We subtract off the variance obtained in this way Irom the 
second moment measured in the simulations, with the results 
shown in Figure 7c. This method ol evolving the shot noise 
tends to overcorrect lor the effects ol discreteness lor the 
earlier epochs and does not lead to a consistent second mo- 
ment on large scales at later times. This can be understood 
as a breakdown ol the approximation made in equation (14). 

In conclusion, the discreteness corrections discussed 
above cannot be applied consistently to the moments ol the 
dark matter particles in the A'-body simulations. However, it 
is important to try to disentangle the effects ol discreteness 
and nonlinear evolution in order to be able to determine the 
scale on which the perturbation theory results break down. 
The approach we follow here is to reduce the problem by 
using simulations with larger numbers ol particles, so that 
the discreteness effects are transferred to scales smaller than 
those ol interest. One can then use the moments as measured 
without applying any corrections. 



4 PERTURBATION THEORY (PT) 
4.1 Linear PT and ^2 

As the initial conditions in our simulations correspond to a 
Gaussian field all the higher order correlation lunctions are 
zero.^ In linear perturbation theory (PT) the Gaussianity 



^ As densities must be positive, this is only possible in the limit 
of small variance, ^2 1. In our particular case, the initial 

arrangements have non-zero ^ j but these are very small compared 
with the ones induced by gravity. 
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Figure 8. The volume averaged 2-point correlation ^2 fo^" 
ferent expansion factors in our largest simulations as a function 
of the comoving radius of the sphere. The continuous lines show 
the theoretical initial conditions, at t = to? and the corresponding 
linear theory prediction scaled to match selected expansion fac- 
tors. The long-dashed and short-dashed lines show the shot-noise 
contribution arising from a grid and from a random distribution 
respectively. 



is preserved and we still have |j = for J > 2. Thus, in 
linear PT the evolution of the density field is completely 
specified by the 2-point correlation, or its power spectrum, 
which grow with the scale factor a: 



(15) 



where x is the comoving scale and to is some initial time. 
The scale factor follows the standard evolution equation: 



3 '' 



(16) 



with k = for the spatially flat models (Q = 1) considered 
in this paper. 

The linear growth in equation (15) is illustrated in Fig- 
ure 8, which shows the values of {2 for different expansion 
factors in our largest CDM simulation (ensemble B). The 
continuous lines show the initial conditions obtained from a 
numerical integration of the initial power spectrum (equa- 
tion 10), using equation (7) in the limit ki = and k2 = 00. 
Each line is normalized using the expansions factors a ob- 
tained from the simulation. 

Figure 8 shows the very good agreement of the theo- 
retical initial conditions with the 1st output time (bottom 
curve) for large scales. For small scales, R < 2/t~^Mpc, 
there is a significant deviation caused by the shot-noise. As 
the initial displacements from the grid are small we can ap- 
proximate the total variance {2 in the 1st output time by the 
direct sum of the fluctuations arising from the grid pattern 

^2 , with the fluctuations imposed by the initial displace- 
J (J 

ments {2 > i-s. equation (14). The long-dashed line in Figure 
8 shows the values of {2 measured for the initial arrange- 
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ment of particles in a grid. After taking into account this 
contribution to the second moment there is good agreement 
between the linear PT prediction and the initial conditions 
for all scales. 

Figure 8 also shows a good overall match between the 
linear PT prediction, equation (15), and different output 
times for R > 5/t~^Mpc. Note that the linear PT pre- 
dictions are normalized with the values of a correspond- 
ing to each output time and therefore both the shape and 
amplitude of the predictions are fixed. The shot-noise con- 
tribution in the estimation of {2 becomes less important 
as we approach the final output time because the intrinsic 
fluctuations grow in amplitude. In these simulations, with 
Npar = 10^, even Poisson shot-noise, shown as a short- 
dashed line in Figure 8, is considerably smaller than the 
intrinsic signal in the last output time. Thus, the deviations 
from the linear theory prediction that can be seen in the last 
output time correspond to intrinsic non-linear evolution in 
the mass distribution rather than to shot-noise. 

The most noticeable discrepancy due to nonlinearities in 
Figure 8 appears at scales R < 5 Mpc where the N-body 
results give larger fluctuations than the linear prediction. A 
smaller discrepancy occurs between R = 8 - 20 h'^ Mpc 
where the N-body result is lower in amplitude than the lin- 
ear prediction (c/the power spectrum measured for this en- 
semble in Paper 1, where we noted that the nonlinearities 
caused a transfer of power from large scales to small scales). 
Although this is a small effect it is quite significant, given the 
size of the dispersion over 10 realisations of the initial con- 
ditions; in Figure 8 we display 2a errorbars. In this regime, 
i.e. i? = 8 — 20 Mpc, neither the shot-noise nor the finite 
size of the simulations make a significant contribution to {2 
(see section §3). 

We can directly check this by comparing the results 

from simulations of different sizes. This is shown in Figure 9 

which compares the last output time in the large simulation 

(i.e. top in Figure 8) with the corresponding output time in 

Ensemble A, i.e. 180 Mpc box with 64'^ particles (see also 

Figure 4). The la errorbars are obtained from the variance of 

10 realisations of the random phases. We also plot the result 

—(2) 

for the second-order perturbation theory {2 obtained from 
a numerical integration, i.e. equation (7), of the second-order 
power spectrum presented in Paper 1. 

Figure 9 shows that there is a very good agreement 
between the two ensembles''' up to i? ~ 10 Mpc. The 
small differences between the two ensembles at larger scales 
are caused by the finite volume of the box. To show this we 
estimate {2 by truncating the linear P{k) at ki = 2Tr/ Lb, 
which corresponds to the largest scale in each ensemble (with 
Lb = 300 or Lb = 180 Mpc). This mimics the lack of 
power at larger scales and thus roughly accounts for the 
finite volume of the simulations. We plot the ratio of this 
last estimate to the full linear prediction with fci = in 
Figure 9, shown by the long-dashed lines. The biggest effect 
occurs for the smaller simulation. Ensemble A. The relative 
change in {2 induced by the finite volume of the simulations 
at i? > 10 Mpc seems to account for the discrepancies 
between the two ensembles. 

+ The corresponding results for Ensemble D, with 32 particles, 
are almost identical to the ones in Ensemble A, see Figure (4). 
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Figure 9. Ratios of the averaged 2-point correlation ^2 to the lin- 
ear PT prediction ^2 ^- Squares correspond to the mean of 10 re- 
alizations of the last output time in Ensemble B, i.e. 300 Mpc 
box and 10^ particles, whilst open symbols correspond to Ensem- 
ble A, i.e. 180 Mpc box and 64'^ particles, starting from a grid 
(triangles) or from a glass (circles). The solid line shows the ratio 
of the second-order to the linear PT. The dashed lines show the 
effect of the finite size of the simulation box. 



Thus, we believe that the discrepancies between the 
simulations and linear theory correspond to intrinsic non- 
linearities that are significant at least up to = 30 Mpc, 

though smaller than 10% for R > ^ h'^ Mpc. The second- 
—(2) 

order PT result, {2 , dashed line in Figure 9, gives slightly 
better agreement but only in a very restricted range of 
scales, i.e. from i? > 15 Mpc. For scales smaller than 
i? = 10 Mpc the second-order result is even worse than 
the linear one. This is not specially surprising (as there is 
a priori no reason why the 2nd order term should make the 
perturbation series converge) and agrees with the results 
found in Paper 1 (see its Fig. 11). 



4.2 Second-order PT and ^3 

To find non-zero theoretical predictions for higher order cor- 
relations, { J, one has to extend the perturbation analysis be- 
yond linear theory. To do this, we expand the perturbation 
equations in powers of the density contrast 8, 

8{X,t) = 8^^\x,t) + 8^^\x,t) + 8^'^\x,t)..., (17) 

where is the linear solution, and = e'(«(^))^ is the 
second-order solution, obtained by using the linear solu- 
tion in the source terms, and so on. For Gaussian initial 
fluctuations, the 3-point function vanishes to linear order, 
= 0, and the lowest order contribution is 
^^(1)^(1)^(2) y Since oc we have that: 
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Figure 10. The averaged 3-point correlation ^3 for different out- 
put times in our simulations (ensemble B) as a function of the 
comoving radius of the sphere. The continuous lines show the 
scaling prediction of equation ( 19) using the second output time. 



where we have used the fact that S^^^ is Gaussian dis- 
tributed. That is, £,:^{t) is given in terms of the square of 



the linear 2-point correlation {2 > 
as a function of time: 

|3(x,t) = 0*13(1:, to). 



and thus it scales with a 



(19) 



We have checked these predictions against the moments 
measured from our simulations. Figure 10 shows {3 for differ- 
ent output times in the largest simulations (Ensemble B), as 
a function of comoving radius R. The continuous line shows 
the scaling prediction, equation ( 19) with the second output 
time sls t = to and using the expansions factors a obtained 
from the simulation § . The most noticeable discrepancies m 
Figure 10 comes at scales R < 5 Mpc where the N-body 
results give larger amplitudes than the PT prediction. The 
agreement is quite good for R > 7 Mpc and it is in- 
teresting to note that we do not find here a significant dis- 
agreement at large scales, in contrast to what we found on 
comparing {2 to linear PT. 

The hierarchical behaviour of {3 in PT can also be ex- 
pressed as 



—PT 

«3 



53 il 



(20) 



where the value of ^3 is set by gravitational instability alone. 
So far we have just checked that the scaling relations work 
for different output times. We now want to check that the 
hierarchical amplitude ^3 estimated from the simulations 



S The first output time corresponds to the Zel'dovich approxi- 
mation (ZA) which, although it has the same scaling as PT, has 
a different proportionality constant, see equation (22) below. 



agree with the ones predicted using PT. Bernardeau (1994a) 
has estimated the hierarchical amplitude, ^3, for {3 in the 
top-hat spherical window case: 



34 

y 

dlo: 



71 



+ 71 
I2 



dlogR' 



(21) 



We have compared this analytical result against the numer- 
ical integration by Frieman & Gaztafiaga (1994) and find an 
excellent agreement. The value of ^3 is not constant, as the 
slope of the linear 7i) changes with scale. Nevertheless, 
71 , and therefore ^3 , is a slowly varying function of R be- 
cause the CDM spectrum is close to a scale-free model. For a 
scale-free, power-law spectrum P{k) oc A:", ^3 is a constant 
with 71 = — (ra -|- 3) (Juszkiewicz et al. 1993). On the other 
hand, for a purely unsmoothed field, i? = 0, W(kR) = 1, 
the normalized skewness is 5*3(0) = 34/7, independent of the 
power spectrum (Peebles 1980). 

—PT 

Note that in comparing the PT predictions {3 from 
equations (20)-(21) with the fully evolved N-body results 
there is the ambiguity of whether to use the linear or the 
non-linear (N-body) ^2- H the the hierarchical expression 
{3 = 53^2 holds beyond PT, one should use the non-linear 
{2 in equation (20). In this case, ^3 in equation (21) gives 
the exact prediction. There is also the further ambiguity of 
whether the linear or non-linear value of 71 should be used. 
This last point makes little difference in our case because, 
at large scales, the shape of {2 changes only slightly as it 
evolves (cf Figure 8). 

In Figure 11, we plot the prediction of ^3 from equation 
( 21), shown by the long-dashed line, using the the shape of 
{2 in the initial conditions to estimate 71 . This is compared 
with the values of ^3 estimated directly from {3 in each 
output time, i.e. ^3 = ^3/^2 with {2 given by either: 

• a) the linear theory value of {2 scaled to the correspond- 
ing output time (continuous lines in Figure 11), or 

• b) the non-linear value of {2 obtained directly from the 
simulation (symbols with errorbars in Figure 11). 

The first output time, labeled ag = 0.24, does not match 
the PT result at all. This is because it corresponds to the 
Zel'dovich approximation which yields a different prediction 
for the value of ^3 : 



53^ 



4-1-71 



(22) 



(see Bernardeau & Kofman 1994). Thus the simulation is 
evolved from the grid to an intermediate state, given by 
the Zel'dovich approximation, with a skewness (and higher 
order correlations) smaller than the corresponding gravita- 
tional ones. The prediction for the Zel'dovich approximation 
is plotted in Figure 11 as short-dashed lines and agrees quite 
well with the first output time (the disagreement at small 
scales comes from the shot-noise contribution). As the simu- 
lation evolves the value of ^3 moves closer to the PT result, 
shown by the long-dashed line. 

Nevertheless, in the last output time, erg = 1, the second 
order PT prediction does not quite agree with the N-body 
result scaled to the linear {2 [case a) above and continuous 
line in the figures]. The continuous line is above the second 
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Figure 11. Amplitudes Ss = ^'^^ different output times, 

iabefed by dg , in our iarge simuiations (Ensembie B) as a function 
of the comoving radius R. The continuous line shows the values of 
Ss estimated using the scaled the linear ^2 while symbols use the 
non- linear ^2 • The short-dashed and long-dashed lines correspond 
to the Zel'dovich and non-linear PT approximation. 



order PT prediction for R < 6 Mpc and below it for 
larger scales. The deviation for large scales is small but sig- 
nificant (given the errorbars ) between R = 8 - 30 Mpc. 

This discrepancy is not very surprising given that we 
have already identified small but significant non-linearities 
in ^2 for ^ similar range of scales. When we normalize the 
values of {3 with the non-linear {2 [case b) above and tri- 
angles in Figure 11], we find a much better agreement for 
i? > 6 Mpc. Using the non-linear {2 corresponds to a 
higher order in the perturbation expansion, but it is not 
clear a priori whether all such higher order terms are in- 
cluded on doing this substitution. The good agreement be- 
tween the predictions and the triangles at large scales in 
Figure 11 seems to indicate that all the relevant terms are 
indeed included with this prescription and that the hierar- 
chical relations go beyond the first contribution in PT. 

The increasing difference between the continuous line 
and the triangles for the different output times in Figure 13 
reflects that non-linearities in {2 increase with time. 

Figure 11 shows a transition from the ZA to the PT 
result as the simulations evolve, i.e. as ag increases. It is not 
clear to start with whether this is a real tendency or just 
an artifact caused by transients from the initial Zel'dovich 
displacements. If the later is true then it would be neces- 
sary to start the simulation earlier, i.e. with a lower value 
of as, in order to obtain an accurate estimation of {3 at an 
earlier time, e.g. at ag = 0.40. However, if the tendency of 
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Figure 12. Same as in Figure 11 with an earlier start, (Tg = 0.10. 



53 to increase with time is real evolution one might expect 
a further increase of ^3 for later times, i.e. for erg > 1. 

We can check this by starting our simulations earlier. 
We have run a new set of 10 large simulations starting from 
erg = 0.10 which corresponds to a expansion factor 2.4 times 
earlier. The results are shown in Figure 12. The first output 
time now gives noisier results because the fluctuations are 
now smaller and more difficult to measure with the counts in 
cells technique (see Appendix for a discussion of the choice 
of mass assignment schemes for smooth particle distribu- 
tions). The second output time, ag = 0.40, shows a better 
agreement with the PT results in Figure 12 than in Figure 
11. Furthermore, the later times are not affected by the ear- 
lier start, as expected if the evolutionary trend of ^3 is an 
artifact as explained above. 

We have further checked this tendency by evolving the 
simulations to a later epoch. For practical reasons we do it 
for the Lb = 180 Mpc and N = 64^ simulations (En- 
semble A). The new results for ^3 are shown in Figure 13. 
The first output time now corresponds to ag ~ 0.16. The 
latest output time, ag = 1.25, shows no change in ^3 (trian- 
gles) compared with erg = 1 (either in this same figure or the 
ones in the larger simulations, figures 11-12). Next we do a 
new set of 10 simulations with an earlier start, ag = 0.10, 
and another set with a later start, ag = 0.26. The results 
for the early start are shown in Figure 14. As expected, the 
ag = 0.40 output give slightly better agreement than be- 
fore while the later times give identical results. The results 
for the late start, ag = 0.26, reproduce very well those in 
Figure 11. All this indicates that to be saved from the ZA 
transient, the simulations should be started a expansion fac- 
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Figure 13. Same as in Figure 11 but for Ensemble A. 



Figure 14. Same as in Figure 13 but starting from (Tg = 0.10. 



tor at least ~ 3 times smaller than the first output time we 
are interested in. 

We have also done a set of simulations starting from a 
glass distribution instead of a grid (see Appendix), with 64'^ 
particles. While there are some differences at small scales in 
the earlier outputs between the glass and the grid simula- 
tions, reflecting differences in the shot-noise, the results for 
the later output times are almost identical (see Figures 16). 



4.3 Higer-order PT and 

In general, the perturbation expansion equation (17) gives: 



2( J-l) 



^j{x,to). 



(23) 



We have checked these predictions against the moments 
measured from our simulations for different output times. 
We find good agreement, within the errors, for all orders 
J = 3 — 10. The results are similar to the ones in Figure 10 
but the errors are quite large. 

The hierarchical scaling also predicts: 



-FT 



(24) 



with characteristic amplitudes Sj set by gravitational insta- 
bility alone (see Peebles 1980, Fry 1984, Goroff et al. 1986, 
Bernardeau 1992, 1994b). 

Amplitudes Sj are estimated from the non-linear corre- 
lations measured in the simulations (as in case b) above). 
The results in the largest simulations are shown for the last 
output time, ag = 1, and for J = 3 — 10 in Figure 15. 

Bernardeau (1994b) has estimated using PT all high 
order hierarchical amplitudes, Sj, in the top-hat window 



case in terms of the logarithmic derivatives jj of {2 • 
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and so on. The numerical estimation of jj becomes very 
uncertain for large values of J but their contribution to Sj 
does not seem very important. Thus we have find that up to 
J = 7 setting jj = for J > 2 yields similar predictions for 
Sj. Here, there is also the ambiguity of whether these values 
of J J have to be estimated from the the linear or non-linear 
value of ^2- Again, we have checked that this makes little 
difference in our case because, at large scales, the shape of {2 
changes only slightly as it evolves (cf Figure 8). In practice, 
it is better to estimate jj from the linear { as any small 
fluctuation is enhanced for the higher orders. 

In Figures 16 we plot these predictions and compare 
them with the values obtained from the simulations: Sj = 
?j/?2 J = 3 — 7. Here we only show up to J = 7 

because the errors for J > 7 are quite large at the scales of 
interest here, R > 7 Mpc. 

We show the outputs time for which ag = 0.40, 1.0 and 
1.25 in the 64'^ simulations for two sets (of 10 different real- 
izations each) of initial particle arrangements; one set on a 
grid (triangles - ensemble A), the other one on a glass (circles 
- ensemble F). The errors (from the dispersion in each en- 
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Figure 16. The hierarchical amphtudes Sj = / *he Lj^ = 180 Mpc, = 64'^ simulations for different output times: 

(Tg = 0.40,1.0,1.25. The continuous line shows the PT approximation. Triangles correspond to the values averaged over 10 different 
realizations of the initial displacements from a grid (the errorbars are the dispersion from the mean). Circles correspond to the averaged 
over 10 different realizations of the displacements from a glass. 



semble) are quite similar for the two cases, and only the ones 
from the grid are shown. For clarity, only points with rea- 
sonable errorbars are plotted. The different initial arrange- 
ments give almost identical results. There is a very good 
agreement within the errorbars with PT for scales larger 
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Figure 15. The hierarchical amplitudes Sj = Cj/C2 *he 
Lb = 300/1-1 Mpc, N = 100^ simulations for (jg = 1.0. 



than R ~ 8h~^ Mpc. In Figure 16 we also plot the values 
for different output times showing that the agreement at 
large scales is not affected by evolution, although there is a 
significant increase of Sj at smaller scales. A similar trend 
is found for earlier times but with slightly larger errors. We 
conclude that, at least during the last three expansion fac- 
tors, Sj are constant at large scales, R ^ 7 Mpc, and are 
not affected by the evolution of clustering, just as predicted 
by PT theory. 



5 DISCUSSION 

When numerical simulations are used to model gravitational 
clustering, it is essential to average over an ensemble of sim- 
ulations (Efstathiou et al. 1985) in order to find the sta- 
tistical significant of the results. We show here that for all 
order statistics the (cosmic) variance from different mem- 
bers of the ensemble is quite significant at all scales, and 
not only at large scales, even when averaged over very large 
volumes (see Figures 5-6). In particular, the variance in S3 
for a box of Lb ~ lOO/j-i Mpc can give an overall error 
in the amplitude S3 of 50%, at a 1-sigma level, for any 
one particular realization. This effect has to be taken into 
account on interpreting the clustering of galaxy catalogs, 
where the cosmic variance is probably even larger than in 
the CDM model (as there is more power on large scales in 
the galaxy distribution than in the CDM model, e.g. Mad- 
dox et al. 1990). Gaztafiaga (1994) found that comparable 
Lb ~ 100 Mpc sub-samples from the APM catalog give 
a large scatter in the amplitude of S3 . This indicates that 
the discrepancies between the CfA or SSRS and the APM 
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are not significant but caused by the sampling fluctuations 
expected in the volume traced by the CfA/SSRS. 

On small scales, shot noise due to the discreteness of 
the particles dominates the signal from the clustering. The 
Poisson shot noise correction is not valid when applied for all 
the mass points in the simulation, unless the initial particle 
arrangement in the simulation was random. As the simula- 
tion evolves, so does the shot noise. We found no consistent 
way to correct for the shot noise. The best approach is to 
use simulations with the largest number of particles possi- 
ble, thus reducing the shot noise amplitude and moving the 
scale at which these effects become important well below the 
particular scales of interest. For the largest simulations used 
in this paper (Lb > 180 /t"^ Mpc and N > 64^), the initial 
particle arrangement does not seem to affect the estimated 
correlations at all on scales larger than the mean inter- 
particle separation; this is despite the fact that visually (see 
Figure A4) the appearance of the voids in the simulations 
started from a grid seems significantly different to that of 
the voids in the simulations started from a glass. 

We found evidence for nonlinear effects in the evolution 
of the variance of the density field, with a transfer of variance 
from large to small scales, in agreement with the result found 
for the power spectrum of fluctuations in Paper 1. 

We confirm previous comparisons of PT and simulations 
(e.g. Juskiewicz et al. 1993, Bernardeau 1994a) and extend 
them to higher orders. However, we find better agreement 
between the simulations and the perturbation theory results 
if we use the nonlinear variance {2 measured from the sim- 
ulations to normalise fj, rather than using the linear vari- 
ance, which strictly speaking is the first non-zero PT result. 
This indicates that, at large scales, the hierarchical relations 
S, J — Sj^2 hold for the fully evolved distribution and the 
values of Sj estimated in PT correspond to the exact ones. 
We have also followed the evolution of in time and find 
that during at least the last ~ 3 expansion factors (Az ~ 2) 
the amplitudes Sj remain unchanged, despite the fact that 
the evolve by large factors, ~ 10"'""'. Although we have 
tested PT for a particular CDM model with Q = 1, the per- 
turbation results for Sj are effectively independent of Q or 
A (Bouchet et al. 1992, Bernardeau 1994a). We have also 
studied models with different initial P{k) and find similar 
agreements (work to be presented elsewhere). 

These are important results because it means that one 
can directly predict Sj from {2 measured in galaxy cata- 
logues, without assuming any particular cosmological (Gaus- 
sian) model (i.e. initial ^2; ^ or Hubble constant), and com- 
pare them with the values of Sj estimated from in the 
same catalogue. 

We show this comparison for the APM catalogue in Fig- 
ure 17. The detailed prediction of Sj involves a knowledge of 
the shape of the variance of matter fluctuations. We assume 
here that this shape is similar to the shape of the variance 
in the galaxy fluctuations, i.e. biasing does not affect the 
shape significantly but might affect the overall amplitude.^ 
Other possibilities, such as non-linear and non-local bias- 
ing, have been studied by Gaztafiaga & Frieman (1994). The 
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Figure 17. Hierarchical amplitudes Sj for J = 3 — 5 from the 
APM (Gaztahaga 1994) compared with the predictions in non- 
linear perturbation theory in the case of no biasing. The long 
dashed lines shows the prediction of the standard CDM model. 
The short dashed curves show Sj obtained when the variance 
is computed from the power spectrum measured for the APM 
Galaxy Survey by Baugh &z Efstathiou (1994a). 



shape of the variance, ^2; the APM is obtained by inte- 
grating the three-dimensional P(k) measured by Baugh & 
Efstathiou (1994a) from the same APM sample. We also 
show the values Sj obtained directly from the higher order 
correlations in the APM (Gaztafiaga 1994). These are esti- 
mated from the angular amplitudes sj{9) with conservative 
errors coming from several sources: the dispersion between 
four different zones of the APM, the merging corrections, 
uncertainties from the selection function and uncertainties 
in the inversion factors (note that part of the errors are 
therefore systematic). There is a reasonable agreement for 
all orders, J = 3 — 5, at large scales R ^ 7 Mpc. This in- 
dicates that clustering in APM is hierarchical and very simi- 
lar to the clustering that emerges from gravitational growth 
of small (initially Gaussian) fluctuations, as predicted in PT 
and found in the N-body simulations. 
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APPENDIX Al: INITIAL PARTICLE 
ARRANGEMENT AND DISCRETENESS 

In this Appendix, we describe how ensembles C, D, E listed 
in Table 1 were run from different initial arrangements of 
particles and compare both visually and statistically via a 
power spectrum analysis, the growth of density fluctuations 
in these simulations. 

The most common way of setting up the initial density per- 
turbations in A'-body simulations is to move the particles 
from a cubic grid arrangement, with the displacements cal- 
culated using the Zeldovich approximation (Efstathiou et 
al. 1985, Zeldovich 1970). This is called the 'quiet start' be- 
cause the grid contains no power on scales other than the 
Nyquist frequency of the particles. The theoretical power 
spectrum of density fluctuations can be represented down 
to smaller length scales than could be achieved using parti- 
cles displaced from random initial positions, especially when 
small numbers of particles are used (c/Figure 4 of Efstathiou 
et al. 1985). The random initial pattern of particles has Pois- 
son shot noise present on all scales. When the simulations 
started from a grid are examined visually however, the un- 
derdense regions that would correspond to voids in the real 
universe are actually populated by mass points that have 
small displacements from their initial grid positions. Simula- 
tions started from a random particle distribution do not suf- 
fer from this problem, but have a different power spectrum of 
density fluctuations on small scales because the white noise 
amplitude dominates the theoretical input power spectrum 
on these scales. 

An alternative to these schemes is to perturb an initial ar- 
rangement of particles that is glass-like, in which all the 
particles try to avoid one another without having a regular 
structure (S. White 1993; private communication). Such a 
distribution may be obtained by running a simulation with 
the signs of the particle velocity updates at each timestep 
reversed, giving effectively a 'negative' or repulsive gravita- 
tional force between the particles. 

We ran a 32'^ particle simulation on a 64'^ grid, starting from 
random positions with zero velocities, reversing the signs of 
the velocity updates at each timestep. A slice one tenth of 
the simulation box thick is shown in projection in Figure Al. 
On scales greater than the mean interparticle separation, the 
particle distribution is subrandom and 'glass-like'; we shall 
henceforth refer to this pattern as a glass distribution. 
Figure A2 shows the power spectrum of the simulation at the 
expansion factors given by the key. We computed the power 
spectrum by tabulating the density of the mass points on 
a 64'^ grid using the cloarf-m-cell charge assignment scheme 
(Hockney & Eastwood 1981) and then taking a Fast Fourier 
Transform. Due to the subrandom nature of the particle 
distribution, it is necessary to use a higher order mass as- 
signment scheme than nearest gridpoint. The cloud-in-cell 
scheme is not accurate on scales around the Nyquist fre- 
quency (kNyquist = whcrc L IS the length of the 
box) of the particle grid. The power spectrum oscillates on 
large scales with expansion factor. This can be explained 
by the particles reaching an arrangement with small density 
fluctuations on large scales at a given time whilst still pos- 
sessing a velocity, so that they overshoot the maximally self- 
avoiding distribution. The dotted line in Figure A2 shows 
the slope of a minimal power spectrum with P{k} oc A;* (see 
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Figure Al. A slice from a 'negative' gravity simulation seen in 
projection. The thickness of the slice is one tenth of the length of 
the simulation box. 



Figure A2. The power spectrum of the simulation run with neg- 
ative gravity, at the expansion factors shown by the key. 



Peebles 1980 §28 for a discussion of how this form for the 
power spectrum arises due to the discreteness of the mass 
particles). Ideally, one would expect the power spectrum in 
the negative gravity simulation to tend to this minimal form. 
However, for wavenumbers smaller than the Nyquist fre- 
quency of the particle grid, shown in Figure A2 to power 
spectrum has a ra = 2 power law index. This could be due 
to the aliasing of power from smaller scales. 
In order to compare different unperturbed arrangements of 
particles as initial positions in A'-body simulations, we ran 
the ensembles listed C-E in Table 1. Figure A3 shows the 
power spectrum of density fluctuations averaged over the 
ten simulations of each ensemble, at three selected output 
times, showing the initial perturbations set up on the parti- 
cle distribution, the evolution of the particle distribution at 
an intermediate timestep and for the final output time. The 
expansion factor of the simulation at each output is shown 
by the side of each curve. The linear theory power spectrum 
for the standard CDM model is shown by the dotted line 
at each value of the expansion factor. The theoretical power 
spectrum is well represented almost up to the Nyquist fre- 
quency of the particle grid in both the simulations started 
from unperturbed positions of a regular grid and a glass. 
The limit on the representation of the theoretical spectrum 
arises because we are in essence using a Particle-Mesh or 
PM code to set up the density fluctuations. The theoretical 
spectrum is Fourier transformed to obtain the corresponding 
gravitational potential, then the particles are moved accord- 
ing to the forces derived from this potential. The dynamic 
range in length of the simulation changes as the particles 
begin to cluster, because they no longer sample the density 
in the simulation cube uniformly. 

For a random arrangement of particles, the theoretical power 
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Figure A3. The power spectra averaged over the simulations of 
ensembles C, D and E, at three values of the expansion factor. 



spectrum is truncated by a white noise spectrum when the 
amplitude of the of the theoretical power spectrum falls 
below the amplitude of the shot noise of the particles, 
Pshot = i/Npar- lu Order to represent the theoretical power 
spectrum more accurately in a simulation started from ran- 
dom initial positions, the initial fluctuations would have to 
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be given a larger amplitude, which means that the simu- 
lations is effectively started at a later time. To impose a 
power spectrum on a random distribution of particles with 
the desired amplitude, it is necessary to subtract the 1/Npar 
shot noise arising from the unperturbed particle arrange- 
ment from the theoretical spectrum. In Section 3.5, we ac- 
tually use an ensemble of simulations started from random 
positions for which the input spectrum did not have this 
correction, instead removing the Poisson shot noise from 
the measured second moment. As the simulations evolve, 
we see the usual transfer of power from large scales to small 
scales (c/ Paper 1). The different growth of the ensemble 
power spectra for high k values reflects the accuracy with 
which the theoretical input spectrum could be represented 
on these scales and the differing amounts of noise arising 
from discreteness on scales around the mean particle sepa- 
ration. 

A quantitative comparison of slices from the simulations is 
given in Figure A4. We show a tenth of the simulation box 
in projection at expansion factors that match those used 
in Figure A3. The simulations were set up with the same 
random phases, so similar structures will form at the same 
location within each slice. Any differences in the visual ap- 
pearance of the particle distribution will be due to the ac- 
curacy with which the theoretical power spectrum could be 
imposed upon the initial distribution of particles. The ap- 
pearance of the voids in the simulations started from a grid 
is very different from that in the simulations started from 
a glass, even though this does not show up in the power 
spectrum. 
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Figure A4. The evolution of the particle distribution in simulations taken from each ensemble with the same initial random phases. 
Each panel shows a slice one tenth of the thickness of the simulation box seen in projection at the expansion factor indicated. 



